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^ ■ We develop and demonstrate two numerical methods for solving the class 

pq . of open cavity problems which involve a curved, cylindrically symmetric con- 

O \ ducting mirror facing a planar dielectric stack. Such dome-shaped cavities 



are useful due to their tight focusing of light onto the flat surface. The first 



O ' method uses the Bessel wave basis. From this method evolves a two-basis 



method, which ultimately uses a multipole basis. Each method is developed 
for both the scalar field and the electromagnetic vector field and explicit 
O ! "end user" formulas are given. All of these methods characterize the ar- 

^ \ bitrary dielectric stack mirror entirely by its 2 x 2 transfer matrices for s- 

r-| ' and p-polarization. We explain both theoretical and practical limitations to 

CLc our method. Non-trivial demonstrations are gi ven, including one of a stack- 

induced effect (the mixing of near-degenerate Laguerre-Gaussian modes) that 
k> ! may persist arbitrarily far into the paraxial limit. Cavities as large as 50 A 

H \ are treated, far exceeding any vectorial solutions previously reported. 
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1. Introduction 

There is currently considerable interest in the nature of electromagnetic (vector) modes 
both for free space propagation [1, 2, -3, 4] and in cavity resonators [-3]. In particular, 
recent advances in fabrication technology have given rise to optical cavities which cannot 
be modeled by effectively two-dimensional, scalar or pseudo- vectorial wave equations [6]. 
The resulting modes may exhibit non-paraxial structure and nontrivial polarization, but 
this added complexity also gives rise to desirable effects; an example from free-space 
optics is the observation of enhanced focusing for radially polarized beams [4]. This 
effect is shown here to arise in cavities as well, among a rich variety of other modes that 
depend on the three-dimensional geometry. The main goal of this paper is to present a 
set of numerical techniques adapted to a realistic cavity design as described below. 

Much work involving optical cavity resonators utilizes mirrors that are composed of 
thin layers of dielectric material. These dielectric stack mirrors offer both high refiectivity 
and a low ratio of loss to transmission, which are desirable in many applications. The 
simplest model of such cavities treats the mirrors as perfect conductors (^tangential = 
0, if normal = 0). lu applications involving paraxial modes and highly refiective mirrors, 
it is often acceptable to use this treatment. The mature theory of Gaussian modes (c.f. 
Siegman [7]) is applicable for this class of cavity resonator. When an application requires 
going beyond the paraxial approximation to describe the optical modes of interest, the 
problem becomes significantly more involved and modeling dielectric stack mirrors as 
conducting mirrors may become a poor approximation. Also, one is often interested in 
the field inside the dielectric stack and for this reason must include the stack structure 
into the problem. 

In this paper we present a group of improved methods for resonators with a cylin- 
drically symmetric curved mirror facing a planar mirror. The paraxial condition is not 
necessary for these methods; tightly focused modes can be studied. Furthermore, the 
true vector electromagnetic field is used, rather than a scalar field approximation. The 
planar mirror is treated as infinite and is characterized by its polarization-dependent 
reflection functions rs(^inc) and Tp^Oi^c) for plane waves with incident angle ^inc- Our 
model encompasses both cavities for which the planar mirror is an arbitrary dielectric 
stack, and cavities for which the planar mirror is a simple mirror (conducting with 
rs/p(6'inc) = —1 or "free" with rs/p^Oi^c) = +!)• The opposing curved mirror is always 
treated as a conductor. It should be noted that, for most modes which are highly fo- 



cused at the planar mirror, (modes which are likely to be of interest in applications), 
this limitation can be expected to cause little error because the local wave fronts at the 
curved mirror are mostly perpendicular to its surface. (Of course, for applications in 
which the curved mirror is indeed conducting, our model is very well suited.) On the 
other hand, the correct treatment of the planar mirror can be a great improvement over 
the simplest model. Applications with both dielectric and conducting curved mirrors 
have been and are currently being used experimentally [5, 8, 9]. 

The methods described here belong to a class of methods which we refer to as "basis 
expansion methods". In basis expansion methods, a complete, orthogonal basis (such 
as the basis of electromagnetic plane waves) is chosen. Each basis function itself obeys 
Maxwell's equations. The equations that determine the correct value of the basis coef- 
ficients are boundary equations, resulting from matching appropriate fields at dielectric 
interfaces, setting appropriate field components to zero at conductor-dielectric interfaces, 
and setting certain fields to be zero at the origin or infinity. In the usual application of 
this method, each homogenous dielectric region is allocated its own set of basis coeffi- 
cients. Our methods use a single set of basis coefficients; the matching between dielectric 
layers is handled by the 2x2 transfer matrices of the stack. 

As dielectric stacks have nonzero transmission, optical cavities with this type of mirror 
are necessarily open, or lossy. The methods described here deal with the openness due to 
the stack and the solutions are quasimodes, with discrete, isolated complex wavenumbers 
which denote both the optimal driving frequency and the resonance width^. While 
the dielectric mirror is partially responsible for the openness of our model system, the 
openness is not primarily responsible for mode pattern changes resulting from replacing 
a dielectric stack mirror with a simple mirror. The phase shifts of plane waves refiected 
off a dielectric stack can vary with incident angle, and it is this variation which can 
cause significant changes in the modes, even though refiectivities may be greater than 
0.99. Generally speaking, the deviation of \rs/p{0[ac)\ from 1 is not as important as the 
deviation of arg(r^/p(6'inc)) from, say, arg(rs/p(0)). 

We develop two general methods, the two-basis method and the Bessel wave method. 
The scalar field versions of both methods are also developed and are discussed first, acting 
as pedagogical stepping stones to the vector field versions. The Bessel wave method uses 
the Bessel wave basis which is the cylindrically symmetric version of the plane wave basis. 
This method is described in Section 3. The two-basis method ultimately uses the vector 
or scalar multipole basis. The multipole basis has an advantage in that it is the eigenbasis 
of a conducting hemisphere, the "canonical" dome-shaped cavity. The unusual aspect of 
the two-basis method is the intermediate use of the Bessel wave basis. The two-basis 
method is developed in Section 4. We have implemented both methods and have used 
them as numerical checks against each other. Various demonstrations and comparisons 



^For many modes, there is also loss due to lateral escape from the sides of the cavity. While our 
model intrinsically incorporates the openness due to lateral escape in the calculation of the fields 
(by simply not closing the curved mirror surface, or extending its edge into the dielectric stack), 
this loss is not included in the calculated resonance width or quality factor, Q. Because a single set 
of basis vectors is used to describe the field in the half-plane above the planar mirror, this entire 
half-plane is the "cavity" as far as the calculation of resonance width is concerned. 



are given in Section 5. Our implementations of all methods are programmed in C++, use 
the GSL, LAPACK, SLATEC, and PGPlot numerical libraries, and run on a Macintosh 
G4 with OS X. Limitations of our model and methods are discussed in Appendix A. 
Appendix B discusses plotting modes that are associated with linear polarization and 
Appendix C specifies dielectric stacks that are used in Section 5. 

2. Overview of the Model and Notation 



\^y&^^ — 




Figure 1: The cavity model. 



A diagram of the model is shown in Figure 1. The conducting surface is indicated 
by the heavy line. The annular portion of this surface extending horizontally from the 
dome edge will be referred to as the "hat brim". The dome is cylindrically symmetric 
with maximum height z = R and edge height z = Zg. The shape of the dome is arbitrary, 
but in our demonstrations the dome will be a part of an origin-centered sphere of radius 
Rs = R unless otherwise specified. The region surrounding the curved mirror will be 
referred to as layer 0. The dielectric interface between layer and layer 1 has height 
z = zi. The last layer of the dielectric stack is layer A^ and the exit layer is called layer 
X. The depiction of the stack layers in the figure suggests a design in which the stack 
consists of some layers of experimental interest (perhaps containing quantum wells, dots, 
or other structures [5, 8]) at the top of the stack where the field intensity is high, and a 
highly reflective periodic structure below. 

At the heart of the procedure to solve for the quasimodes is an overdetermined, com- 
plex linear system of equations. Ay = b. The column vector y is made up of the 
coefficients of eigenmodes in some basis B. The field in layer is given by expansion in 
B using these coefficients. For a given wavenumber, fc, a solution vector y = yhest can 



be found so that \A{k) y — 6p is minimized with respect to y. Dips in the graph of the 
residual quantity, A^ = |^(A;)2/best ~ ^1, versus k signify the locations of the isolated 
eigenvalues of k (theoretically A^ should become at the eigenvalues). The solution 
vector yhest{k) at one of these eigenvalues describes a quasimode. The system of equa- 
tions is made up of three parts (as shown below): Ml equations, M2 equations and an 
arbitrary amplitude or "seed" equation. 



Ay 



Ml 

M2 
[s. eqn.] 



y 



■ff 



(1) 



Henceforth Ml refers to the planar mirror and M2 to the curved mirror. 

The Ml boundary condition for a plane wave basis is expressed simply in terms of 
the 2x2 stack transfer matrices Ts{6i^f.) and Tp(^im,), as suggested by the Ml region 
(enclosed by the dashed line) in Fig. 1. In the scalar and vector multipole bases, a sort 
of conversion to plane waves is required as an intermediate step. The dashed k vectors 
in the figure (incoming from the bottom of the stack) represent plane waves that are 
given zero amplitude, in order to define a quasimode problem rather than a scattering 
problem. The plane waves denoted by the solid k vectors have nonzero amplitude. 

The M2 boundary condition is implemented as follows. A number of locations on the 
curved mirror are chosen (the "X" marks in Fig. 1). The width of the hat brim is Wb as 
shown. An "infinitesimal" hat brim {wh ^ A) is introduced to give the dome a diffractive 
edge. An "infinite" hat brim can be introduced theoretically and can make the model 
more easily understandable in certain respects. More about the model in relation to the 
hat brim is discussed in Appendix A. 2. The M2 equations are the equations in basis B 
setting the appropriate fields at these locations to zero. For a problem not possessing 
cylindrical symmetry, these locations would be points. The simplification due to this 
symmetry, however, allows these locations to be entire rings about the z axis, specified 
by a single parameter such as the p coordinate. Finally, the seed equation sets some 
combination of basis coefficients equal to one and is the only equation with a nonzero 
value on the right hand side (6). 

The cylindrical symmetry of the boundary conditions allows one to always find solu- 
tions which have a dependence of exp(«?7i0), where m is an integer. This in turn leads 
to a dimensionally reduced version of the plane wave basis called the Bessel wave basis, 
in which each basis function is a superposition of all the plane waves with the same 
wavevector polar angle, 6k. The weight function of the superposition is proportional 
to exp(zm0fc). We will refer to the non-reduced basis as the "simple plane wave basis". 
The unadorned phrase "plane wave basis" (PWB; same abbreviation for plural) will re- 
fer henceforth to either or both of the Bessel wave and simple plane wave bases. When 
using the scalar or vector multipole basis (MB; same abbreviation for plural), cylindrical 
symmetry allows the problem to be solved separately for each quantum number m of 



interest^. The dimensional reduction in this case amounts to the removal of a summation 
over m in the basis expansion. 

The refractive index in region (layer) q is denoted hq. Layers are also denoted with an 
upper subscript in parenthesis: E^'^' means the electric vector field in layer q. Sometimes 
"fs" is used as a value of g, meaning "in free space" (e.g. rifg), whether or not any of the 
layers in the model we are considering actually are free space. We note here that fifg = 1 
only in "cavity type I" discussed below. 

The symbol k, where not in a super/subscript and not bold nor having any su- 
per/subscripts, always refers to what may be called the "wavenumber in free space", 
although it will have an imaginary part if Ml is a dielectric mirror. An imaginary part 
in wavenumber, refractive index, and/or frequency is often introduced (as it is in this 
problem) to turn open cavity problems into eigenvalue problems. The definition of k is 
as follows. Define —k"^ as the constant of separation used to separate space and time 
equations from the wave equation for layer q: 

V^X(.. ^ Y-^. (2) 

Here X*-^-* may be a vector or scalar field. In a few steps, the selection of a global 
monochromatic time dependence exp(— ztut) reveals that the ratio kq/fiq is independent 
of q. Then k is defined as /c = /cfg, so that kq = riqk where riq = hq/hfy. In the model, 
the index ratios riq are assumed to be real. The single plane wave solution to (2) has the 
form X''^^ = C*^''''e*'" "^ -x^-uvt ^]^g],g (j{q) jg g^ constant vector or scalar and the complex 
wave vector k'^'''^ is given by fc^''^ = kqfl)^ = kuqflj^ with H)^ being the unit direction 
vector of the plane wave, specified by 6j^ and 0^. The generally complex frequency is 
given by u; = ckq/hq. At a refractive interface, the angle 9^ changes as given by Snell's 
law. 

To understand the meaning of a complex fc, it is helpful to realize that the spatial 
dependence of the quasimodes are identical^ in the following two physical cavities: 

I: nfs = l,ng = nq, 

II : fifs = T, fig = TrZg, T G C. (3) 

(The shape and size of each dielectric and conducting region are the same for cavities 
I and II.) Cavity I is composed of of conductors and zero-gain regions of real refractive 
index. Cavity II is constructed by taking cavity I and multiplying the refractive index of 



^The modes with low \m\ are likely to be of practical interest since they have the simplest transverse 
polarization structure. The |m| = 1 family of vector eigenmodes is exceptional because the pro- 
portionality of Ep and i?<^ to exp(im(/)) means that modes for |r7T,| ^ 1 have no average transverse 
electric field, even instantaneously. (It is straightforward to show that (R.eEx),t> = (R.eEy)^ = if 
\m\ ^ 1.) Thus a uniformly polarized, focused beam centered on the cavity axis can not couple to 
cavity modes with |m| 7^ 1 ! 

^There is the minor difference for the magnetic field first noted in Eqn. (12). Since the discrepancy 
is a constant factor multiplying H , however, this difference is not necessarily part of the spatial 
dependence. 



each region, including free space, by an arbitrary complex number T. The congruence 
of the spatial quasimodes follows from separating the variables in (2). The values of fc, 
kq, and Uq for a given quasimode are the same in cavities I and II. The frequency in 
cavity II is Uu = Wi/T. If we henceforth consider only the specific cavity II for which T 
= (1 — ig) where g is tuned to be the ratio {—luik/Rek) (for a given quasimode), we see 
that ujii = cRe/c = Recui. While in cavity I the quasimode decays in time, in cavity II 
the quasimode is a steady state because the gain exactly offsets the loss. Either of the 
two cavity types may be imagined to be the case in our treatment. The only difference 
is the existence of the decay factor exp(clm(/c)t) for cavity I. (The inequality ImA; < 
always turns up for an eigenvalue problem with conducting and/or dielectric interface 
boundary conditions.) We note that the relation of k to the free space wavelength (always 
real) of a plane wave is A; = (1 — ig)2n/X. The quality factor, Q, of the quasimode is 
Refc/(2|ImA;|) = 1/(2^). 

We note that Snell's law, n^sin^^'^ = Uq^ismOj^ , is independent of whether we 
have a cavity of type I or II because the quantity (1 — %g), if present, divides out. One 
of the limitations of our method is the omission of evanescent waves in layer and 
in layers where Uq > Uq (see Appendix A.l). Snell's law may cause 9)^ to become 
complex for layers with index ratios Uq less than Uq. In this case sin 6'^'' > 1 and 
cos^i'^^ = 2Sgn(cos^f )[sin2^[:^^ - 1]V2. 

In most cases the symbols ip, E, and H stand for complex-valued fields. The time 
dependence is exp{—iujt) and it is usually suppressed. Physical fields are obtained by 
multiplying by the time dependence and then taking the real part. 

Throughout this paper, the common functions denoted by Yim, Pr, Pi, Jn, jh ^ind n/ 
are defined as they are in the book by Jackson [10]. 

In the implementation, c and the related constants eo, /io, and Zq are all unity, and 
they will usually be dropped in our treatment. We also assume non-magnetic materials 

so that jjLq = jjLq = 1. 

3. Plane Wave Bases and the Bessel Wave Method 

Although this major section describes the Bessel wave method, much of what is discussed 
here is applicable to the two-basis method with little alteration. The discussion in Section 
4 is greatly shortened due to this overlap of concepts and procedures. 

3.1. The Field Expansion in the Simple Plane Wave Bases 

3.1.1. Scalar basis 

A single scalar plane wave in layer q has the form ij) = Cexp (ik^''^ • x — icot). For a 
general monochromatic field, k and u are fixed and the field can be expressed (due to 
the completeness of the scalar PWB) uniquely (due to the orthogonality of the scalar 
PWB) as a sum over plane waves in different directions. In our treatment however, we 
omit plane waves in layer q which would only exist as evanescent waves when refracted 



into layer 0. The expansion for the field in layer q is 



^(i)[x) = r d0fc r dOf sm{ef)^^^e'^'''--. (4) 





Here the basis expansion coefficients are the 'ij)^ (continuous coefficients in the integral, 
and discrete coefficients in implementation). The above expansion effectively propagates 
each plane wave existing in the cavity down (whether forward or backward) into the stack 
layers, and adds up all of their contributions. In order to express the ip^ in terms of 
ipl.\ it is first necessary to separate the coefficients with k^ > from those with k^ < 
and write the above expansion as 

V;(«)= f\(f>, r daf sin(4°^) 
Jo Jo 

X f^^'^^e*''"''-"' + ^pfe'^^"'^-'^) . (5) 



The u and d refer to the plane waves going upward or downward, e.i. tpu is the expansion 
coefficient ipj^ for ki > {or kl > 0, since sgn{kz ) = sgn{kz )) and ipf takes the 
place oiijjj^ for ki' < 0. The wavevector fc*^^^ in cylindrical coordinates is (fcp , k^f , kz ), 
for which the following relationships hold: 

A;('?) = /tgsin^^^ = knosine'^^ = knosma^^'> = kf, 
k-z — kqCosOj^' = knqSgn{cos9l ) coso;^ 

= A;r2gSgn(cos0^ )W 1 — I — sina^ 



riq 



M 



This leads to 



(indep. of q). (6) 



p27v /'7r/2 

Jo Jo 

x^^We^^^+^i-^^e-^^^), (7) 



where 



^p = pkriQ sin(a[°^) cos(0 - 0^), 

Lp^ = zkuq COS a^ . (8) 

From standard theory regarding plane waves and layered media[ll], one can calculate 
the 2x2 complex transfer matrix, TJ , that obeys the following equation 



Defining the column sums +Ps = T^u + ^s 22 ^^^ +1^ = T, 1 1 + T;2i allows us to 
write the expansion of the scalar wave in the layers as 



xe'^''(+/3i'?)^W++7i'^)V^f). (10) 

The reason for the notation with the subscripts "+" and "s" will become apparent in the 
vector discussion. The "s" refers to s-polarization. 

The variables in the simple PWB for scalar fields are the complex ^i and ^^ (the 
superscript will often be dropped). Next we consider the simple PWB for vector fields. 

3.1.2. Vector basis 

We assume that for our purposes a general monochromatic electromagnetic field can be 
expressed uniquely as a sum of vector (electromagnetic) plane waves. For every given 
frequency and wavevector direction Qk there are two orthogonally polarized plane waves 
(as opposed to a single plane wave in the scalar case). Instead of a single coefficient ijjk for 
each spatial direction we need two, Sk and Pk, which we can define as follows. Sk is the 
amplitude of the vector plane wave propagating in direction fi^ which has its electric 
field polarized in the x-y plane {E^ = 0). Thus, this plane wave is an "s-wave" with 
regard to the planar mirror. Pk is the amplitude of the "p-wave", the vector plane wave 
in direction Qk which has its electric field polarized in the plane of incidence {E^ = 0). 
The coefficients Sk and Pk will be separated into Su, Sd, Pu, and Pd- 

To specify the polarization of the fields we will use unit vectors denoted by e. The 
unit vector e^"^^ denotes the direction of the electric field associated with the plane wave 
with wavevector k^''^ and s-polarization. We take the direction of the unit vectors to be: 

= xsiiKpk -ycos(f)k, 



e 



(9) _ fl(9)o„„/'„„c,fl(9) 
p,k 



6>^''^sgn(cos^^^ 



= pfcSgn(cos0^ ) cos6'^ — 2sgn(cos6'|. ') sin6^^' . (11) 

In this phase convention (used by Yehfll]), the projections of the e^l vectors for the 
incident and refiected waves onto the x-y plane are equal. The other common phase 
convention has these projections being in opposite directions. 

The entire electric and magnetic field can be broken up into two parts: E^'^^ = Es + 
Ep and -H"*^^^ = Hs' + Hp where Hs is the field with magnetic s-polarization 
{Hz = 0) and Hp is the field with magnetic p-polarization {H^ = 0). We can now 
write down the most compact expansion of the vector field. 



4^) = / dnfPt^^e^^le^'^'''- 



10 



/f (^) = n, I dnf^Si'^el^nicos ef)e^''-. (12) 

The factors of fiq in the H equations come from the physical relation oi H to E for a 
plane wave. Note that hq is different for cavity types I and II (as given in Eqn. (3)). 
Separating up and down coefficients yields 

Ei^^= Td^^ei^l r^'dafsin(4°))e^-^ 
Jo Jo 

l-2n i'Tt/2 

El,"^ = / d0fc / da[°^ sin(a[°)) e^'^" 
Jo Jo 

+ z sin(4^^) f-P^'^h'^^ + P^'^^e-^^A] . (13) 

These expressions explicitly use coordinate vectors only where necessary due to a de- 
pendence of the efc vectors on the sign of cos 6'[, . The expressions for i?*^*^^ are omitted 
for brevity. 

To relate 5^% and P^V^ to S^^^ and P^^ we can use the transfer matrices: T^ for 
s-polarized light and Tp for p-polarized light. The transfer matrix used for the scalar 
field in Eqn. (9) is the same matrix we will use here for s-polarization. These matrices 
perform the following transformations 

Pfe- 







We define 



i/Ji-^) ^ Tit ± t/'^) 



s,225 



±fs — -^3,21 =■= -^5,11; 

n{q) ^ rp{q) , rpiq) 
±fjp — ■lp^i2 =L -^p,22) 

±7i^)=T(t±r(J. (15) 

Note +/3s and +7s are defined as before. The P and 7 quantities are functions of z 
and Zi and not of p or 0. They are functions of k and a^ but not of 0^. 



11 



Eii) 



dal,°^ sinfal"^^-*^" 



Now the field expansions become 

r"27r rn/2 

'0 



p,cos(4'^))(+/?(«)pW++7i'')p^ 



(<?) p(0) 



+ z sm( 



in(«i:^^)f-^i^^pi°^--7S^^^r 



i/(^) 



n„ 



27r 



d0.^S 



7r/2 



dai°^ sin(a[°))e*'^'' 



x(-/5^^)Pi°)--7i^)pr), 



2n 



n/2 



//!«) = n„ / d0fc / da[.°^ sinf«l°h e*'^" 







X 






(16) 



The variables in the simple PWB for vector fields are the complex 5*1 , 5*^ , P^ , and 
PJ (the superscript will often be dropped). 



3.2. The Field Expansion in the Bessel Wave Bases 

3.2.1. Scalar basis 

We have already assumed a time dependence of exp(— «ci;t). As mentioned in the 
Overview, a cylindrically symmetric set of boundary conditions allows us to assume 
an azimuthal dependence of exp(?m0) with m being an integer. Consider the expansion 
(4). We wish to find the conditions on ipj^ which cause the entire dependence of il>{x) 
on (j) to be exp(zm0). 



The general Fourier series expansion of ijjj^ is 



^?H^r\0.)=E/^' 



{9)|'/^(o)^ f,^ri4>k 



(17) 



We can then write (4) as 



^(^) = y: / dot 



W^in^flW^.^^^fWr/gW 



sm 



^^^.fiM 



27r 



5fc e 



ipkno sm{e]^')cos{if>-<j>k)pin(l>k 



(18) 
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where 

(Pz = zkuq cos 9 j^ . (19) 

The last integral is of the solved form 



'0 

where J„ denotes the regular Bessel function of order n {n can be negative). This yields 

^(5) = 27r V(z)"e^"<^ r defsm{ef) 

xe^^^J„(pA;nosin4°V^'n4°^ (21) 



In order to have only exp(zm0) dependence on </>, the integral on the right hand side 



must be zero for n ^ m. Because fn {9\. ) cannot be a function of z or p, the only way to 



have this for all z and p is to pick /„ = for n 7^ m. Thus '^k\0f\ 0fc) = /m Vr)e™"^^ 
At this point we define the symbol il)^ to mean the coefficient fm ■ The cylindrically 
symmetric expansion is 

^("^ = e r de^^h\n{ef^)e^^^J^{pkn,s\r.ef)i,^\ef^), (22) 

where 

i = 27T{i)"'e'""t'. (23) 

This is an expansion in scalar Bessel waves, defined to be 

^ exp [izkuq cos O]!' ) J^ (pA;no sin 6*^ ^ ) , (24) 

with {wk \^k )} being the set of coefficients. Each Bessel wave is a set of simple plane 
waves with fixed polar angle but having the full range (0 to 27r) of azimuthal angles, 
(pk- The weight factors of the plane waves are proportional to exp(«m0fc). The final 
cylindrically symmetric scalar expansion with up and down separated is 

-i/;('?)(a;) = ^ / daf^ sin(af^) J„(pA;no sinaf^) 
Jo 

x(+/5i^Vf ++7i^Vf). (25) 

The il>u and ^^ are the (complex) variables in the Bessel wave method for scalar fields; 
they make up the solution vector y in (1). 
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3.2.2. Vector basis 

For cylindrically symmetric boundary conditions, the dependence of Ep, E^, E^, Hp, 
Hz, and H^ can be taken (for a single mode) to be exp(zm0). 

Consider doing the 0^ integrations in (12) or (16). The unit vectors ^^/ f, and Pk 
depend on 0^. The z components do not depend on (pk so we will look at these first. 
There is no contribution to the z component of the electric field from Eg nor is there 
any contribution to the z component of the magnetic field from Hg . We define 

fP^"?) (x) = /f («) . z = H^p"^ ■ z. (26) 

Requiring that these quantities have an exp(zm0) dependence produces results similar 
to the scalar case. Defining 

oil) ^im(t>k — oiq) 
'-'k ^ — '^k ^ 

P^'^^e'""^'' = P^''\ (27) 

and using (16), we have the final, useful expansions for ^P^'^^ and ^P^''^: 

^p{q) = ^ r dftW sin(c,(0)) sin(«^'')) 
Jo 

fP^'^= -^n, r^'daf sin(«i°))sin(«i:')) 
Jo 

X J^pknosinaf^) [+Pi'^Sl^^ + +ii'^Sf^) . (28) 

To deal with the transverse part of the electromagnetic field it is helpful to use quan- 
tities related to circular polarization. We define 

±^(«) = ±tEi'i^ ■ X - E^f^ ■ y 

= e^"^{±tEi''^ ■p-E'f^ -(j)), 
IS'^'^'^ = ±iH'f'^ ■ X - H^"^ ■ y 

= e±^^±zifi^).p-/fi'').0), 
±P(«) = ^i") • X ± z^i^) • y 

= e±^^4^).p±z4'').0), 

Hp{q) _ jj{q) . ^ ^ ^jj{q) . ^ 

:/f(^).p±zif('').c^). (29) 



^±^4 



Inverting (29) yields 



—I 
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^W . ^ = Z!( ^pWg-*'^ - „p('?)e*'^), (30) 

with the magnetic field quantities having similar relations. Now we use (12) and (11) 
with (29). The resulting electric field quantities are 

^S^i) = r d^f sin(ef )e'^^ 
Jo 

(\(h,Q^'^'t'k g«pfc"0 sin{e[, ) cos{4i-(j)k ) Q 


JO 

r-27r 



Jo 

p(i) = f d^f sin(^f)cos(e[''^)sgn(cos^f )e*^^ 
Jo 

i'2tt 



It is the exp(±z0fc) factors in the integrands here that motivated the definitions of ±S/P 
(29). We see that the substitution of (27) into (31) results in 0^ integrals of the form 
(20). Performing this substitution, doing the integrals, and separating the up and down 
parts gives the final expansions: 



±5^''^ = ^± / daf' sm{af')Jm±i{pknosma)^') 

/.7r/2 

.^pW = ^_j_ / da^^' sm{a^^')Jm±i{pknosma]^') 
Jo 

xcos(4^^)(+/3(^)pi°)++7i^)pr), 

/.7r/2 

±S^'^^ = $,±hg / daf^ sm{a^^')Jm±i{pknosma^^') 

/.7r/2 

^p(9) = ^j_n^ / daf^ sm{af')Jrn±iipknosmaf') 
Jo 
,(.))^__^(,)^(o)+_^(.)5(o))^ (32) 



X cosl^a^. 



where 



^± = 27r(z)'"±^e*('"±^)'^. (33) 



At this point one can quickly verify, using (30) and the above equations, that Ep, E^ 
Hp, and iJ^ do indeed have a 0-dependence of exp(zm0). 
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The Su , S^j^ , Pu , and PJ are the (complex) variables in the Bessel wave method 
and make up the solution vector y in (1). They are essentially coefficients of electro- 
magnetic Bessel waves, although we need not explicitly combine (32), (30), (28), and 
(26) to obtain an explicit expression for the E and H vector Bessel waves as we did for 
the scalar case (24). 

3.3. The Linear System of Equations for the PWB 

Until now the PWB coefficients have been treated as continuous, when in practice they 
must be chosen discrete. Let us keep in mind this discrete nature in the following 
subsections. We denote by A^dirs the number of directions a\. we choose. Thus there 
are 2A'dirs coefficient variables for a scalar problem and 4A^dirs coefficient variables for a 
vector problem. The distribution of the a}, on [0, 7r/2] need not be uniform, and the 
effect of distribution choice will be briefly mentioned later. 

3.3.1. The planar mirror (Ml) boundary equations 

The Ml equations (planar mirror boundary condition equations) in the PWB are very 
simple. In fact, because of this simplicity, the Ml equations in the MB are basically a 
transformation to and from the PWB with the Ml equations for the PWB sandwiched 
between. The reflection of a plane wave off of a layered potential is a well known 
problem. For the purpose of determining the field in layer 0, the entire dielectric stack 
is characterized by the complex rg and r^ coefficients acting at the first surface of the 
stack. For the scalar case with the layer 0-layer 1 interface at ^i = 0, the boundary 
condition is just ^i = ^^(a^ )^^ where rs{a) is the stack reflection function. Since 
Bessel waves are linear superpositions of many plane waves with the same 9l parameter, 
this equations is true for Bessel waves: 

^i°^=r.(«r)^r- (34) 

For a conducting mirror, set r^ = — 1 and for a free mirror set r^ = 1. 

If the interface is at a general height 2i, then the same r^ function is used at this 
surface yielding 

or 

V^f - rs{af^)i^P = 0, (35) 

where 

Vp(«r) = V,(ar)e-^2-^"— i^ (36) 
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(The equation is given for both s- and p-polarization since we will shortly be using the 
rp quantities.) The quantities rs/p(a^ ) are independent of z and zi. When M2 is a 
dielectric stack mirror, the fs/p{a) are determined by T^J~ ^^'^^ according to 

n/pia) = -T2Jt2,2- (37) 

For the vector case the equation for the s-polarized plane waves is 

E^^-efl = Uaf)Ef-e% (38) 

where E}^ is the total electric field of the two plane waves going in the direction specified 
by afc, 0fc, and u or d. Shifting to our current notation and to Bessel waves, the equation 
becomes 

S^^-Uaf)sT = ^. (39) 

For p-polarization there is an arbitrary conventional sign. In our phase convention 
(chosen in (11)) the equation is 

^i°^ • e^Sl = fMf)ET ■ eZ, (40) 

or 

^i°^-r>r)^f =0. (41) 

The vector e^ ^ ,^ is e^ '- with k being k forced into the up/down version of itself. 

For a simple mirror, set rp = r^ = {—1 for conducting, +1 for free} and use (36) 
instead of (37) to determine fg/p. Sometimes we will use cosak as the explicit argument 
to fg/p or rs/p instead of a^. 

Equation (35) or Eqns. (39) and (41), given for each discrete a}. , form the Ml 
equations. 

3.3.2. The curved mirror (M2) boundary equations 

As mentioned in Section 2, the M2 equations come from setting the appropriate field 
components equal to zero at some number of locations on the curved mirror and the hat 
brim. If we chose individual points on the two dimensional surface, the 0-dependence 
factor, exp(2m0), could be divided away. Thus picking locations with the same p and z 
but different yields identical boundary equations. Therefore we simply set = (and 
t = 0) and pick equations by incrementing a single parameter (such as p) on the one 
dimensional curve given by the intersection of the conducting mirror and the x-z plane. 
The number of locations, Nm2 loc, determines the number of M2 boundary equations. 
All of the locations are taken to lie in layer 0. 

We obtain the M2 equations by, in effect, doing the a^ integrals in (25) or in (28) and 
(32). Before making the integrals discrete, the distribution of representative directions 
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must be chosen. If the interval Aa[. between successive directions is not constant, the 
integral over a^ must be transformed to an integral over a new variable, x, where Ax 
is constant. Such a transformation will generate a new integration factor. At this point 
the integral is turned into a sum according to: f^ -^ ^ ■ , dx ^ {b — a)/N. Choosing 

the direction distribution to be uniform in a^. requires no change in integration factor 
and yields a "summation factor" of 7r/(2A'dirs)- 

Using (25) and (2.3), the M2 equations for the scalar problem are 



27r I ^^ j Yl ['^n^{pMo sin a'j^. 



^dirs/ j^^ 



xsin(4°))f+/3f|_z^., + +7f|_^, 



■j y 1 -r- i.- lz=z^ r u,j 



0. (42) 



An equation is added to the linear system for each chosen location specified by (p*, z*). 
All phase factors have been divided out of (42) but the scale factor vr^/A^dirs has been kept 
for representative weighting. Of course there is also an effective weight produced by the 
distribution of the evaluation locations on the conducting mirror. In our implementation, 
we choose equal steps of 6 to cover the dome and equal steps of p to cover the hat brim 
(see Fig. 1). 

For the vector problem there are three equations associated with each location: E^ = 
0, E\\ = 0, and if^ = 0. (Here the subscript "||" refers to the direction that is both 
tangential to the M2 surface and perpendicular to (f).) From (.30) the E^ = equation 
is 

— ( +5(°) + _5(°) + ^( +P(°) - _P(°))) = 0. (43) 

The expansions for ±S/P in terms of the unknowns S/ P^L in equation (32) must now 
be used, along with the identical integral-to-sum conversion used in the scalar problem 
(42). It is probably not beneficial to work out the long form of this boundary equation, 
as its computer implementation can be done with substitutions. 

The E\\ = equation depends on the shape of the mirror. If r] is the angle that the 
outward-oriented surface normal makes with the z axis, then E\\ is given by 

-E|| = Ep cos r] — Ez sin r], (44) 

where E^ = zP^^^ and, using (30), 

Ep = ^ (z( _5(°) - H_5(°)) + +P(°) + „PW) . (45) 

Again equation (32) and the integral-to-sum conversion must be used to obtain the 
explicit row equation. The H± = equation is obtained by doing the same type of 
substitutions with 

H± = Hp sin T] + Hz cos 77. (46) 
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Here H^ = ^P^^^ and 

For locations on the hat brim, r] = 0. 

3.3.3. The seed equation 

All the Ml and M2 equations have no constant term. Thus the best numerical solution 
will be the trivial solution ?/best = 0. To prevent from being a solution, an equation 
with a constant term must be added. One simple type of equation sets a single variable 
equal to 1, for instance Su{j=5) = 1. Another simple type sets the sum of all of the 
coefficients equal to 1. A more complicated type sets the field (or a field component) 
at a certain point in space equal to a constant. No one type of seed equation is always 
best. 

3.3.4. Solution of Ay=b 

As depicted in (1), the matrix A is made up of the left hand sides of the Ml, M2, and 
seed equations. For the scalar case there are 2A^dirs columns and (A^dirs + ^m2 loc + 1) 
rows. For the vector case there are 4Adirs columns and (2A'dirs + 3A'm2 loc + 1) rows. The 
value of A^m2 loc is picked so that A has several times as many rows as columns. A value 
of k is picked and the overdetermined system of equations is "solved" as well as possible 
by a linear least squares method. The best such methods rely on a technique known as 
singular value decomposition [12]. Our implementation relies on the function zgelsd of 
the LAPACK fortran library. To find the eigenvalues of fc, the imaginary part of k is 
set to zero and the real part of k is scanned. As mentioned in the Overview, this results 
in dips in the value of A^. Using Brent's method [12] for minimization, the minimum 
of the dip is found. The real part of k is now fixed and Brent's method is used again 
to find the best imaginary part of k. Then Brent's method may again be used on the 
real part of k. By this alternating method, the complex eigenvalue of k is found, along 
with the eigenmode, y. In practice Brent's method need only be used two to four times 
per scan dip to get an accurate complex k. We usually normalize each row of A to 1 
so that the normalized error per equation in the system can be expected to be around 
A„ = Ar/[|y| X (number of rows in A)^/^]. A„ is one indicator of the accuracy of the 
solution. 

3.3.5. Calculating the field from y 

Once y is calculated for a quasimode, the values of the field in any layer can found by 
using the expansion (25) for the scalar field and equations (-32), (.30), (28), and (26) for 
the vector field. Of course the integrals over a^ must be made discrete as discussed 
previously. Appendix B explains more regarding the plotting of the fields. 
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4. Multipole Bases and the Two-Basis Method 

As mentioned in the Introduction, the two-basis method ultimately uses the scalar or 
vector MB. The MB is the eigenbasis for the closed, conducting hemisphere or sphere. 
Both the vector and scalar multipole bases have known forms which we will use but 
not derive. The basis functions already possess an azimuthal dependence of exp(zm(/)) 
and the dimensional reduction due to cylindrical symmetry is accomplished by picking 
a value for m instead of summing over basis functions with many m. 

The method of stepping along a one-parameter location curve to obtain the M2 equa- 
tions is the same as for the PWB. Explicit formulas in the MB of course will be completely 
different and will be given in this major section. The methods of solution to the linear 
system of equations are the same as for the PWB. The development of the Ml equations 
in the MB, however, requires considerable work. After the system of equations has been 
solved, using the resulting solution vector y to calculate/plot the fields in layers other 
than layer also requires significant work. We use the term "two-basis method" because 
of the role of plane waves in these two calculations. Figure 2 represents the linear system 
of equations of the two primary methods and how they are related. 



JTMm L^J, f IL« r lJ 




Bessel wave 

lTH!(hfjd 



two- basis 
rtielhod 



Figure 2: Diagram for the two primary methods. The closed loops suggest the self- 
consistency or "constructive interference" of the quasimode solutions. Grey 
regions indicate intersection between PWB and MB methods. Size roughly 
indicates the "post-basis-derivation work" required to get the equations. The 
variable coefficients for the vector problem are shown. 
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4.1. The Scalar Multipole Basis 

The scalar MB functions we use are the ipim = ji{knor)Yim{0,(j)) where ji denotes the 
spherical Bessel function of the first kind and F/m is the spherical harmonic function. 
The scalar MB functions, like the scalar PWB functions, satisfy the wave equation. We 
assume the field in layer 0, in a region large enough to encompass the cavity, can be 
expanded uniquely in terms of the scalar MB functions. Using the cylindrical symmetry 
of the cavity to solve the problem separately for each value of m, we expand the field in 
the cavity as 

iP^'^Xx) = J2 ciji{knor)Yim{0, 0). (48) 

l=\m\ 

The expansion coefficients, q, are complex. One should never need to choose /max much 
larger than Re(fc)nor'jnax where rmax is the maximum radial extent of the dome (not the 
hat brim). (Semiclassically, the maximum angular momentum a sphere or hemisphere of 
radius Rg can support (for a given k) is ~ Re{k)noRs, which corresponds to a whispering- 
gallery mode.) 

The scalar MB functions are the exact eigenfunctions of the problem of a hypothetical 
spherical conductor specified by r = .R^ = _R, with eigenvalues given by the zeros of 
ji{knQR). The basis functions for which (/ + m) is odd are the eigenfunctions of the 
closed hemispherical conductor. This is because 

Yim. has parity (-1)'+™ in cos^ (49) 

and thus is zero at 6^ = 7r/2. It can also be shown, using the power series expansion of 
Pi{x) given in Ref. [13], that YimiT^/^, (p) is nonzero if {I + m) is even. 

4.2. The Ml Equations in the Scalar MB 

The expansion of a plane wave in terms of the (monochromatic) scalar MB functions is 
[10]: 

oo I 

e^fc- = 47r^ 5^ {i)%U0.'P)3i{kr)YU0k.<Pk). (50) 

/=0 m=-l 

The inverse of this relation, the expansion of a scalar MB function in terms of monochro- 
matic plane waves, is 

ji{kr)Yim{e, <P) = j A^k {^^Yi^iOk^ <Pk)\ e*'^-". (51) 

This equation is easy to verify by inserting (50) into the right hand side. The use of 
Yim = (— 1)*"^*-^ cind the orthogonality relation for the Yim leads directly to the left 
hand side. The expansion (51), applied to layer (k -^ ko = krio), is the foundation for 
this section and its counterpart for vector fields. 
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^(0) ^ / ,^(0) ^..(0).. ^ c^-^Y^UOfA.) ■ (52) 



Using (48) and (51) yields the scalar field expansion for a given m 

\i 

\=|m| ^ 

The quantity in the curved brackets is ifjk, the simple plane wave coefficient (compare 
to Eqn. (4) with q set to 0). We can now use (35). For k^ > 0, ipu = i^k and ipd = V'fc'i 
where in cylindrical coordinates k = k (fc) = {kp,k(f,,—kz). For kz < 0, ipd = i^k and 
^u = A'- Using (49), one gets ^.^ = E; Q>'^m(^r, 0fc)(-l)'+™(-^)7(4vr). We can now 
derive an equation that expresses (35) and holds for all Jl^ : 



, 1 — rs, / + misevenl , , 

X ....... .(0)^ ,.^ , ,^._.. =0. (53) 



l=\m.\ 

_ sgn(cos 9]^ ) {1 +fs) , / + m is odd J 
At this point there are two ways to proceed. 

4.2.1. Variant 1 

One way to construct the Ml portion of the matrix A is to pick some number A^mi dirs 
of discrete directions, 9], \ and use (53) for each direction (the (pk dependence divides 
out) . Inspection of the equation reveals that it is even in cos 9\, ; thus only polar angles 
in the domain [0, 7r/2] are needed. Using a\. as before to denote this reduced domain, 
the Ml equations in this variant become: 

E-^^>^"^(«r,o) 

l=\m\ 

X ('l - (-l)'+™r,(cos4°^)e-*2.^'="°™'°^°') = 0. (54) 

We note for later comparison that for each a\. we need only evaluate a single associ- 
ated Legendre function (inside the F/m), because the recursive calculation technique for 
computing P™ (x) also computes PJ^^x) for \m\ < I < /max- Since each step of this 
recursive calculation involves a constant number of ffoating point operations, we may 
say that the complexity associated with the P{^ calculations for each a], is 0(/^ax)- 
Generally Nyn dirs ~ /max so the P™ complexity is 0(/max) ^^d the overall number of 
evaluations of fg is 0(/max)- 



4.2.2. Variant 2 

Rather than picking discrete directions to turn (53) into many equations, we could 
project the entire left hand side onto the spherical harmonic basis, {Yii^i}, (that is, 
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integrate (53) against Yjf^, (J7^, ) in f2^ ). Due to the uniqueness of the basis expansion, 
the integral for each [I ,m) pair must be zero. The integrals for m ^ m are zero and 
can be neglected. Thus the number of equations generated is A'^,' = l^^^ — \m\ + 1. We 
^pt /' — / 

^^^ ''max "-max- 

Since (53) is even in cos 6*^ , integrating against Y^f^ gives if / + m is odd, halving 
the number of equations. When / +??7, is even the integration must be done numerically. 
The most analytically simplified version of the coefficients of q in the Ml equation 
corresponding to / (for / + m even) is: 



Ml,', = c f p!"'\x)plr\x) (1 - (-i)'+™f.(x)) dx, 

Jo 



(55) 



where 



47r V^ '^ ^(/+ |m|)! (/ + |m|)!y ^ ' 

This comes from using the properties oiYim under a sign change of ?7i, using the definition 
of Yim: doing the 0fc integral, and halving the domain of the even integral over x = 
cos al . Noting that the number of unknowns in the solution vector y is Ni = /max — 
\m\ + 1, the number of (complex) coefficients that must be calculated is about Nf /2. 

Noting that Ml^;' = (— «)' ~'M1;' , when / + m and / +m are both even, the calculation 
is reduced to about ?>Nf /A real-valued integrations. 

While this variant of the method is in some ways the most elegant, the numerical 
integrals are extremely computationally intensive. Our implementation used an adaptive 
Gaussian quadrature function, gsl_integration_qag of the GSL. (Adaptive algorithms 
choose a different set of evaluation points for each integration and achieve a prescribed 
accuracy; in a sense the integral is done in a continuous rather than a discrete manner.) 
For an adaptive algorithm, the complexity'* associated with the P™ calculation is 0{l^^^, 
where u >1. The number of evaluations of f^ is 0(/^5^)- I^ practice, variant 2 is much 
slower than variant 1 even for cavities as small as R/\ ~ 5 with simple mirrors. Checking 
between the two variants has generally shown very good agreement. 

4.3. The Linear System of Equations in the Scalar MB 

The Ml equations have been given in the previous section. Using (48) yields an M2 
equation 

tmax 

5^Qjz(/cnor,)yi„,(e„0)=0, (57) 

l=\m\ 



^The exponent v is defined so that the number of integration points that must be sampled (for the 
most complicated integrals where Z ~ / ~ ^max) to maintain a constant accuracy goes like ZJ^^x- 
(This definition may not be rigorous if ^ is a function of /max-) Since Pi{x) has 0{l) zeros and the 
P™ for TO > are even more complicated, it is reasonably certain that v >1. 
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for each location (r*,6'*). The discussions from Section 3.3 regarding the M2 equations, 
the seed equation, and the method of solution to Ay = b apply here. The number of 
columns in A is Ni and the number of rows is {Nyn jirs + -^m2 loc + 1) for variant 1 or 
about {Ni/2 + Nm2 loc + 1) for variant 2. 

4.4. Calculating the Field in the Layers with the Scalar MB 

To calculate the complex field anywhere in layer once a quasimode solution y has been 
found, Eqn. (48) can be used directly. To calculate the field in the layers below the 
cavity (g > 0) where the expansion does not apply, the Bessel waves must be used, along 
with the Ts matrix that propagates them. Performing the 0^ integration in (52), using 
(49), and then comparing with (9) and (25) with g = in both of these yields the Bessel 
wave coefficients: 

YUc^fM-ir-. (58) 

I 

Now (25) can be used with g > yielding 

f7r/2 

2 Jo 

X Jmipkno sin a^°^) ^ ci{-iyYir,^{a^^\ 0) 



u 


~v *- 


■T 


/L^ 47r 



^(.)(a,) = i!i_e^-'^ / ' daf sin(af^ 



I 



X (+/?(") + (-l)'+"V7(''))]. (59) 

This is a costly numerical integration to do with an adaptive algorithm. For every 
sampled a^ , the quantities Tg, P™ and Jm must be calculated once. For displaying 
large regions of the field in the stack, it is sufficiently accurate to simply convert the 
solved MB solution vector y into a PWB solution vector (with iVdirs ~ ^max) by means 
of a separate program using (58), and then use the discrete form of (25) to plot. 

4.5. The Vector Multipole Basis 

We expand the electromagnetic field in layer (at least in a finite region surrounding 
the cavity) in the vector multipole basis using spherical Bessel functions of the first kind 
(ji). The multipole basis uses the vector spherical harmonics (VSH; same abbreviation 
for singular), which are given by[14] 



Mi^{x) = -ji{knor)x x VYi^{e,< 



Nimix) = —V X Mi^. (60) 
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The VSH are not defined for / = 0. The nature of the electromagnetic multipole ex- 
pansion is developed in section 9.7 of the book by Jackson[10] using somewhat different 
terminology. The multipole expansion of the electromagnetic fields is 

^max 

'' ^niin 

^max 

/f W (x) =fioJ2 i'^i^im + hiNi^) , (61) 

^ 'rnin 

where cylindrical symmetry has been invoked to remove the sum over m, and /min = 
max(l, \m\). The ai and bi are complex coefficients and there are A'^; = /^nax — ^min + 1 of 
each of them. The ai coefficients correspond to electric multipoles and the hi coefficients 
correspond to magnetic multipoles. The explicit forms of the VSH are 

^ / tTTl 

Mim{x) = —^ji{knor)Yi^{e, 0) 

- ( d 

+ (l)i-ji{knor)—YtmiO, 

^ / 7777 r) \ 

+ '^ I ^^ (rUknor)) ¥^6, 0) . (62) 

\ kn^r sm 6 or J 

4.6. The Ml Equations in the Vector MB 

The goal of the calculation here is to transform the Ml relations, (38) and (40), into 
equations such that each dot product is written in the form ^i{aifa{k) + hfh{k)). The 
resulting (two) equations will be the vector analogues of (53) and will hold for all fi^ , 
leading to two variants of solution method in the same manner as before. First we must 
Fourier expand £J^°^ as we expanded 7/)^°^ to get (52). The only Fourier relation we have 
is (51) which expands the quantity jiYim{x). We must therefore break Mi^ into terms 
containing this quantity. 

Here we can make use of the orbital angular momentum operator L = —i{x x V). 
Using the ladder operators, L± = L^ ± iLy, the VSH Mim = —ijiLYim can be shown to 
be 



Mira = {-tjl) X 



*2 i^tXl^+l + dimyim-l) 



I 



+ :y 2 (~^«m^im+l + di,J'lra-l) + z (mY,„)J , (63) 

where 

dt ^ V(^Tm)(/±m + l). (64) 
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Using (51) multiple times yields 



with 



dnfm^e^'''-, 



(-7V f-?y 



^ . . lm~ y Q^ K ^J'^'lm 



\l 



+ Z- 



A-K 



-niYimi^f), 



.(0)^ 



,(0)^ 



where 

We can avoid the task of expanding Nim into jiXira terms by just using 



V X M,^ / do(o)^fc(o) X M,^e^^'°^-. 



^ /cnn 



Performing the cross product and substituting into (61) using (65-66) yields 



(65) 



(66) 



(67) 



(68) 



EJx) 



Ey{x) 



(0) ikW- 



d^l" e 



^^(0)^.fc(0).. 
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47r 



ai ( -m sin 91 sin 0fcYi„, - - cos O^ ' e,^ ) + 6; ( -e+„ 
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3(0) 



3(0) 



3(0) 



ai [ m sin 6*^ cos 0^1^^ - ^ cos 61^ e^;^ 



'/ -TT^ 



Irn 



B,W=/dni°'e-»-{$:(=l)i[a,( 



'^ ■ i(o)_. ^ .- , l.;^Z)(0) •, 



ai ( - sin 6'^ cos 0^ e,„ + - sin Ol sin 0^ e^ j + hi [mYim 

(69) 

where Yim is understood to mean Yim{'^\)- The quantities inside the curly braces in 
(69) will he denoted hy E^^k), Ey{k), and Ez{k), where the "(0)" superscript on k has 
been dropped. We will also need E{k ) where k (fc) = {kp, k^, —k^) as before. Again 
using (49) we find 



/ 



I 
I 



l+m 



l+m 



l+m 



ai[-m sin 6*^ ^ sin (pkYim - - cos of^ e^^ 
ai[msin9fcos<P,Yi^-\cos9fhl 



bii^ret. 



bi[^ei^ 



ai [ — - sin ef^ cos (pk e;„^ + — - sin 6^^ ^ sin (pk e^^ ) +hi{ niYir^ 



(70) 
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4.6.1. The Ml equations for s-polarization 

We will now define fim and gim by 

\i 
47r 



E{k) • e,,fc = J2 V^ ^"'^'™ + ^'^'™) • (^1) 



z 
Performing the dot product using (11) yields 

-msmef'Yim 



\{^l{^^^ -"^^i^-"^- m-i^m+ie-'"^ 

+ ./{l + m){l + m-l)Yi.,^m-ie"^A • (72) 



The simplification in the last step has been done using recursion relations for the Pl^. 
Performing the dot product for E{k ) yields 

E{k') ■ e,,, = J2 ^(-1)'+" {aigin. - bifim) • (73) 



I 



An 



The Ml relation (from (38)) is E{k) • es,k = rsE{k') • e^^fc for A;^ > and E{k') • e^^k 
fsE{k) • Es^k for kz < 0. These lead to an Ml equation that holds for all J7^ : 



tin ax / 

-I 



E K-i'j I aisigim -h oiS2jim, L -h m even i ^ , . 

A-K \aiS2gim + biSifim, l + moddj ' 



y_ jaiSigim + biS2fim, I + m even\ 

A-TT 

where si = 1 — fs{\ cos6'[, |) and S2 = sgn(cos6'|, ^)(1 + fs{\ cos6'[, ^|)). 



4.6.2. The Ml equations for p-polarization 

The Ml equation for p-polarization is obtained the same way as above. The calculation 
is simplified by temporarily switching phase conventions. We define the unit vector 

Ep^k = sgn(cos6'^ )^p,k (using (11)) so that 

ip^k = X cos 9^^' cos (f>k + y cos 9^^' sin (pk- z sin 9^^' , 
^pk' ~ ~ ^ ^'^^ ^k ^'^^ 't'k ~ y cos 9f, sin (pk — zsm9j,' . (75) 

The Ml relation is now E{k) • Sp^k = {-rp)E{k') • e^,^' for k^ > and E{k') • e^^^i = 
{—fp)E{k) • ip^k for kz < 0, where fp has not changed in any way. Performing the dot 
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products, we find that 

An 



E{k) • Ep^k = ^ -J — {-aifim + hgir 



I 



^(^') • S,fc' = E ^(-l)''-'"(«'/'™ + ^i9im). (76) 



I 



This leads to the Ml equation for p-polarization which holds for all i7^ : 

E {-iY \-a1P2fim + hPigim, I + m eveni ^ ^ .^^. 

47r \-aipifim + biP29im, l + moddj ' 

where pi = 1 + fp{\ cos6'^ |) and j92 = sgn(cos^[. )(1 — fp(| cos6'^ |)). 

4.6.3. Variant 1 

Variant 1 for the vector case proceeds exactly as it does for the scalar case. The Ml 
equations (74) and (77) are seen to be even in cos O^ , allowing us to pick directions 
only in the domain [0, 7r/2]. Thus for each discrete direction a^ we get one equation (in 
the Ml portion of A) for each polarization. No modification of (74) or (77) is needed, 
except that we can take 0fc = and cos6'j. = cosa^. > 0. 

The complexity of calculating the Ml portion of A is of the same order in /max as it is 
for the scalar MB. However, numerous constant factors make the computation time an 
order of magnitude longer, as there are more Pp^ functions to compute and fp must be 
computed in addition to fg. 

4.6.4. Variant 2 

In this minor section, we assume that m > (see Appendix B for plotting negative m 
modes from positive m solutions). We follow the same procedure used in variant 2 for 
the scalar field. Integrating (74) and (77) against Yi'm'i^k ) yields zero if either m ^ m 
or / + m is odd. For m = m and / + m even, the Ml equation associated with each / 

in {\m\, \m\ + 1, . . . , /max} for s-polarization is 

'max ] [ r^ / \ / \ 

E -24^7o (i-(-i)'^'"^"O(^'"^'^'^^ + ^^+"'^(^+"'"^)^"^'('^0^'^('^)^'^ 

^min 

/I 
(^1 + (-1)'+"V,) (^Pr^\x) - (/ + m)(/ - m + l)Pr-\x)^Pp{x) dx = 0, 

(78) 



1=1 
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and the Ml equation for p-polarization is 



'max 1 r /"^ / \ / \ 

J2-2^ a^-^)/ [l-{-iy'-"'fp){^Pr+\x)-{l+m){l-m+l)Pr~\x))Pp{x)dx 

^ ^min 

+ bi (l + (-l)'+'"fp) (P;?+'(a;) + (/ + m)(/ + m - l)P;'!.V'(2;))^r(^) ^^ 



0, 

(79) 



where P" 



for |n| > fi. Note that if m = the equations for / =0 are included. We 
implemented this by calculating an array of each of the following types of integrals: 

J^ p™+ip- dx, J^ P;r-'^r da:, 
J^Pjn+ipjj^rsdx, /,'p--ip-f,dx, 

fl pm+lpm^ J rl pm-lpniif; J™ 

Jo M V P ^■^' Jo /i ^ u ' P '-'•^• 

The integrals in the first row can be calculated ahead of time and stored in a data 
file. The other integrals must be calculated each time the matrix A is calculated. The 
complexity of the calculation is of the same order in /max as for variant 2 in the scalar 
case. These integrals are very numerically intensive and in practice variant 2 takes 
much longer than variant 1. We have found good agreement between the two variants, 
indicating that variant 1 can be used all or most of the time. 



4.7. The Linear System of Equations in the Vector MB 

The M2 equations are found in the usual way, choosing locations (p*,^*) with = 
and constructing the equations E^ = 0, E\\ = Epcosr] — EzSinr] = 0, and H± = 
HpSinr] + H^cosrj = 0. The fields E and H are found by substituting (62) into (61). 
The derivatives in (62) are 



d_ 
dr 



{rji{knor)) = knorji_i{knor) - lji{knor), 



(80) 



and, for m > 0, 
d 



de 



Yin 



(2/ + l)(/-m) !]'/' 
47r(/ + m)! 



X 



(p, 



m+l 



(/ + m)(/-m + l)P™-^)e™^ (81) 



where P/^^ 



0. (For negative m use YJ^-m = (~l)'"^im-) The number of columns in 
A is 2Ni and the number of rows is (2A^mi dirs + 3A^m2 loc + 1) for variant 1 and about 
{Ni + 3Nm2 loc + 1) for variant 2. 
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4.8. Calculating the Field in the Layers with the Vector MB 

To calculate the field in layer 0, the expansion (61) is used directly, with (62), (80), and 
(81) being used to obtain the explicit form. As for the scalar case, there is no direct 
expansion for the layers g > in the MB, and we must use the Bessel wave basis. Using 
(71), (73), and (76), the conversion is seen to be 

'max {_^Y 



Su= ^ — {aigim + kflm) 



n(0) (0) , „ 



s^=j: 



HV, 






{-iy+"^{aigirn-bifi^) 



i"'=-i'\'P^=o 



Pu — 2_^ ~^ {—O-lflm + hgirr^ 



I 



'i:'^="i"^<^fe=o 



P'^ = Y. V^(-l)'+"' («J^™ + bWlr^ 



All 
I 



M) (0) , „ 



(82) 



This is the vector case analogue of Eqn. (58). 

At this point the Bessel wave coefficients Su-, etc., are continuous functions of a^ , 
not discrete as they were in the Bessel wave method. These continuous coefficients are 
substituted into the integrals (32), and finally (30) is used to give the vector fields in the 
layers. As in the scalar case, however, doing the integrals with an adaptive algorithm is 
impractical if one is plotting the vector fields in a sizable region. In this case, it is best 
to make the Bessel wave coefficients discrete and use a simple integration method, or, 
as was suggested in the scalar case, to pick a iVjirs and create a PWB solution vector 
that can then be plotted via a PWB plotting routine. 

5. Demonstrations and Comparisons 

In this section we will demonstrate, but not thoroughly analyze, several results that we 
have obtained using our model and methods. The modes shown here all have A„ < 
0.0002 and do not contain large amounts of high-angle plane waves, which in principle 
can cause problems (see Appendix A.l). The authors expect to publish a separate 
paper focusing on the modes themselves. In this section we will also compare our 
implementations of the two-basis method and the Bessel wave method. Before reading 
this section it may be helpful to look at Appendix B. 

5.1. The "V" Mode: A Stack Effect 

As mentioned several times in the Introduction, the correct treatment of a dielectric 
stack can be essential for getting results that are even qualitatively correct^. One of 



^It is not always essential however, see Note Added in Proof. 



30 



%; 




Figure 3: Contrast-enhanced plot of ReE^ in the y-z plane for a "V" mode {m = l,k = 
8.1583 - 20.000019). M2 is spherical and origin-centered; R = 40; zi = 0.3; 
Ze = 3; Ml = stack II; kg = 8.1600; Ns = 20. In this and other side view plots, 
the field outside the cavity has been set to zero. 

the most remarkable differences that we have observed when switching from simple to 
dielectric Ml mirrors occurs right where someone looking for highly-focused, drivable 
modes would be interested: the widening behavior of the fundamental Gaussian mode 
(the 00 mode). A Gaussian mode of a near-hemispherical cavity will become more and 
more focused (wide at the curved mirror and narrow at the flat mirror) as Zi is decreased 
(from a starting value for which the Gaussian mode is paraxial) . As the mode becomes 
more focused, of course, the paraxial approximation becomes less valid and at some 
point Gaussian theory no longer applies. We have found, using realistic Ali„a;Gaa;As- 
AlAs stack models (stacks I and II described in Appendix C), that the 00 mode splits 
into two parts, so that in the side view a "V" shape is formed. Figure 3 shows a split 
mode and Fig. 4 shows the physical transverse electric field, ReEx = ReE^x + ReEyy, 
for this mode near the focal region. Figure 5 shows the values of Re{kR) as zi is changed 
from 0.3 to about 0.63. Figure 6 shows a zoomed view of the field at the focus of the 
mode at zi = 0.63, where it is qualitatively a 00 mode. If a conducting mirror is used, 
the central cone simply grows wider and wider, but does not split. The V mode is 
predominately s-polarized and appears in the scalar problem as well. Thus it appears 
that the V mode is a result of the non-constancy of arg(rs(a)) for a dielectric stack. 
We note here that following the 00 mode as zi is changed is an imperfect process: it 
is possible that the following procedure may skip over narrow anticrossings that are 
difficult to resolve. However, the character of the mode is maintained through such 
anticrossings. 

The apparent splitting of the central cone/lobe for higher order Gaussian modes has 
also been observed. In the next section we look at some higher order Gaussian modes, 
focusing not on what occurs at the breaking of the paraxial condition, but on what is 
allowed and observed for modes that are very paraxial. 
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Figure 4: An x-y cross section of the mode in Fig. 3 at z = 1.58. Units are in /im 
as usual. For m = +1, the forward time evolution is simply the counter- 
clockwise rotation of the entire vector field. If m = +1 and m = —1 are mixed 
to create a cosine mode, the arrow directions remain fixed and their lengths 
simply oscillate sinusoidally in time so that the mode is associated with x 
polarization. The inset shows the location of the cross section in the view of 
Fig. 3. (The field is plotted only in layer 0.) 

5.2. Persistent Stack-Induced Mixing of Near-Degenerate 
Laguerre-Gaussian Mode Pairs 

As mentioned in the previous section, we do observe the fundamental Gaussian mode 
for paraxial cavities with both dielectric and conducting planar mirrors. The situation 
becomes more interesting when we look at higher order Gaussian modes. It is true 
that, as modes become increasingly paraxial, Gaussian theory must apply. However, 
the way in which it applies allows the design of Ml to play a significant role. Here 
we will demonstrate the ability of a dielectric stack (stack II) to mix near-degenerate 
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Figure 5: Following kR as Zi is varied. The points of the two graphs correspond one-to- 
one, with the parameter zi increasing from left to right in both graphs. The 
peaks in resonance width have not yet been analyzed. 

pairs predicted by Gaussian theory into new, more complicated near-degenerate pairs 
of modes. We must provide considerable background to put these modes into context. 
Our discussion applies to modes with paraxial geometry, modes for which the paraxial 
parameter h = \/{nwo) = wq/zr is much less than 1. (Here wq is the waist radius and 
zj^ is the Rayleigh range, as used in standard Gaussian theory.) In this section, the 
solutions we are referring to are always the +m or — m modes discussed in Appendix B, 
and never mixtures of these, such as the cosine and sine modes discussed in the same 
section. 
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Figure 6: The field in layers 0-X for the mode at zi = 0.63. The inset shows the 
entire mode. The "V" nature of the mode has been lost as it has become 
more paraxial and more like the fundamental Gaussian. KcEt everywhere lies 
nearly in a single direction at any instant (here it is in the x direction). The 
x-z plane is shown here, although the plots of E^ in any plane containing the 
z axis are very similar. Here k = 8.2098 — zO. 0000401. 

5.2.1. Paraxial Theory for Vector Fields 

From Gaussian theory, using the Laguerre-Gaussian (LG) basis, we expect that the 
transverse electric field of any mode in the paraxial limit is expressible in the form: 



Eq 



N 

Eh: 

j=0 



t( 1 



A: 



^^min{j,N-j)(^)- 



(83) 
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Here A^ > is the order of the mode. The Af are complex coefficients and ( ^ ) and 
(ij) are the Jones' vectors for right and left circular polarization, respectively. The 
explicit forms of the normalized ^Gj^^,-j^_-^ functions are given in Ref. [15] as "wj^m" 
with the substitutions n — > N—j and m — > j. The important aspects of the LG^^j"/ ^_--. 
functions are that the ^-dependence is exp(«(2j — N)(j)), and the p-dependence includes 
the factor L'^^Z. 'j^_.J2p'^/w{z)'^) where L^ is a generalized Laguerre function and w{z) is 
the beam radius. In the paraxial limit the vector eigenmodes of the same order become 
degenerate. There are 2(A^ + 1) independent vector eigenmodes present in the expansion 
(83). 
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Table 1: Vector Laguerre-Gaussian modes. 
Independent from the discussion above, Et = Epp + E^(f) can be written as 



Ej 



^-{Ep + zE^)e^^l^^^)+^{Ep-zE^ 



-i4> 



(84) 



oc cxp (im4>) oc cxp (im4>) 

Since the E'p and i?,/, fields have a sole ^-dependence of exp(zm0), comparing with (8-3) 
reveals that, for a given m, at most two terms in (83) are present: the A~ and A'^, 
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terms where j^ = {N + m ± l)/2. Explicitly, 

{E, ± zE^)/2 ^ Aj^LG|S-Un(„.±i,-,.Ti)]/2- (85) 

If the solution (84) has both terms nonzero for almost all x, the constraints < j < A^ 
force A^ to have a value given by iV = \m\ + 1 + 2z/, i/ = 0, 1, 2, . . . . However, if only right 
(left) circular polarization is present in the solution, A^ = |m| — 1 is also allowed, provided 
that m > +1 (m < -1). So, given a (paraxial) numerical solution with its m value, we 
can determine which orders the solutions can belong to. The reverse procedure, picking 
A^ and determining which values of m are allowed and how many independent vector 
solutions are associated with each m, can be done by stepping j in (83) and comparing 
with (84) or (85). The results for the first four orders are summarized in Table 1. The 
LGo({) and LGo(ij) modes are the fundamental Gaussian modes. 

Since the cavity modes are not perfectly paraxial, the 2(A^ + l)-fold degenerate modes 
from the Gaussian theory are broken into A^ + 1 separate degenerate pairs. The truly 
degenerate pairs for m 7^ of course consist of a +m and a — m mode which are related 
by reflection (see Appendix B). The pairings are shown in the last column of Table 1. 

The dashed boxes in Table 1 enclose pairs of modes which may be mixed in a single 
solution for fixed m (Eqn. (84)). Only the mixable m = modes are exactly degenerate 
and may be arbitrarily mixed. For the other mixable pairs, the degree of mixing will 
be fixed by the cavity, in particular by the structure of Ml. We now discuss our results 
regarding the mixing of the two modes with N = 2 and m = +1. 

5.2.2. A demonstration of persistent mixing 

Figure 7 shows the cross sections of two near-degenerate m = 1 modes found for a 
conducting cavity. Here zi = Ze = and M2 is spherical with radius Rg = 70 but is 
centered at z = —59.5 instead of the origin, so that R = 10.5. The paraxial parameter 
h is at least as small as 0.1 (see the side views in the inset plots). Mode A corresponds 
very well to the pure LG]^(|) mode, while mode B corresponds very well to the pure 
LGo(ij) mode. These are mixable modes, as indicated in Table 1, and the conducting 
cavity has "chosen" essentially zero mixing. Note that mode B would have zero overlap 
with an incident fundamental Gaussian beam centered on the z axis, while mode A 
would have nonzero overlap. 

Pictures C and D in Figure 8 are cross sections of near-degenerate m = 1 modes for a 
cavity with stack II. To show that the transverse part of these modes are mixtures of A 
and B, we have added and subtracted the solution vectors j/a and j/b to create the new 
(non-solution) vectors: 

Vc' =yA + vvb, Vy)' =yA- (0.225)77^3, (86) 

and have plotted these fields (C and D ) using k = {k^ + k-Q)/2. (The scaling factor 
?7, here ~ 0.07, is unphysical and related to the effect of the seed equation on overall 
amplitude.) Comparing C and D with C and D shows that we have reconstructed the 
mode cross sections quite well (up to a constant factor). Modes C and D are not pure 
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Figure 7: 8 x 8 /im cross sections of modes of a conducting cavity near maximum 
amplitude {z = 0.25). A;a = 7.89285 and ks = 7.89290. The inset plots show 
KeEx in the x-z plane with horizontal and vertical tick marks every 1 /xm. 

Hermite-Gaussian (HG) modes, but their resemblance to the HG02 and HG20 modes is 
not an accident. Mode conversion formulas in Ref. [15] show that LGq is made of the 
HG02, HGii, and HG20 modes while LG? contains only HG02 and HG2o- Note that both 
modes C and D would couple to a centered fundamental Gaussian beam. 

An interesting property of the mixed modes C and D is that their general forms 
appear to be persistent as h is varied, provided that the paraxial approximation remains 
sufficiently applicable. Furthermore, we do not have to be particularly careful about 
setting Xs to correspond very closely to 27r/Re/c. The modes shown in Figure 9 are not 
nearly as paraxial as the C and D modes. Here ks = 8.1600 which is not close to the 
k values of the modes. Furthermore, stack I was used which is missing the spacer layer 
of stack II. Nevertheless the modes of Fig. 9 bear a strong resemblance to the more 
carefully picked C and D modes. 

5.3. Modes with m 7^ 1 

The m = modes may be circularly polarized, such as the paraxial m = modes of 
Table 1. The more interesting polarization, perhaps, is the m = analogue of x and y 
polarization: radial and azimuthal polarization. The mode shown in Fig. 10 is radially 
polarized, although it is impossible to tell this from the plots because the left and right 
circularly polarized modes are radial (identical to the plot) at t = and azimuthal 
(like a vortex) at ut = 7r/2. The radial or azimuthal polarizations are easier to obtain 
numerically than the right or left polarizations because radially (azimuthally) polarized 
modes have all of the bi (ai) MB coefficients equal to zero and thus can be selected by 
the seed equation (see Appendix B). 



37 





t-1 '.'j ■>■'+"■■''" ■■i."'*-»r''-'— 
' ^ * ■* ' - ' - -■■ ^^■■- 



m> ■ 






•;-'.%v.v^v^i*; -' :^^^- ^ ' :: H ■ ^->-'- -J.^'V v^ 










-■ 






»T- <-'_"_*- T.* ■ j'-'' r ■■ i"- -TT-'' .'T'-'f ■'¥''■'"- " -"-I- 






n^^^n 


T. 




^^1 


. 



Figure 8: Cross sections of true modes (C and D) of a cavity with stack II near maximum 
amplitude {z = 0.05), and constructed modes (C and D ) near their maximum 
amplitude {z = 0.25). The constructed cross sections are nearly identical to 
the true cross section, kc = 7.74685-z0.00005415, kj^ = 7.74680-z0.00005516, 
and the constructed modes are plotted with k^' ,^r = 7.8929, the average of /ca 
and /cb- The inset plots show ReE^ in the x-z plane. The stack parameters 
are N, = 22 and k, = 7.746814. 
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Figure 9: Cross sections of neighboring, poorly paraxial modes near maximum am- 
plitude {z = 1.33). M2 is spherical and origin-centered with R = 10. Ml 
is stack I with A^^ = 20, kg = 8.1600, and Zi = Ze = 1.0. The mode 
on the left has k = 8.51160 — zO. 0002491 and the mode on the right has 
A; = 8.51540 - iO.0003184. 

Figure 11 shows a m = 2 mode that corresponds to the LGi( j) mode from Gaussian 
theory. We also have completed the N = 2 near-degenerate family by finding the LGq ( I ) 
(m = 3) mode for both the cavity used in Fig. 7 {k = 7.8926 ^ kA ^ ks) and the cavity 
used in Fig. 8 (/c = 7.7466 — 20.00005466 ^ kc ^ k^)- The cross sections of this mode for 
the two cavities, which are not shown to conserve space, appear identical, as predicted 
by the absence of a m = +3 near-degenerate partner for this mode (Table 1). 

5.4. Nd\rs and Imax- Comparing the Two Primary Methods 

The authors of this paper first implemented the more complicated two-basis method, 
believing that many more PWB coefficients than MB coefficients would be required to 
expand dome-shaped cavity fields. When the Bessel wave method was implemented as 
a check, it was discovered that the PWB works surprisingly well, with usable values 
of A^dirs being of the same order of magnitude as usable values of /max- The choice 
of seed equation, A^dirs or /max, and other parameters can both affect the depth and 
narrowness of the dips in the graph of A^ versus Rek. This makes it difficult to perform 
a comprehensive and conclusive comparison of the two methods. 

One practical problem with the Bessel wave method is that its performance is some- 
times quite sensitive to the value of A^dirs- In both methods, setting the number of basis 
function too high results in a failure to solve for ^best; the program acts as if Ay = b 
were underdetermined, even though the number of equations is several times the num- 
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Figure 10: Tightly focused, radially polarized m = modes of a conducting cavity. 
The focusing has caused E^ to be greater than the Ep. The spot size of the 
dominant E^ field for such modes can be surprisingly small, as has recently 
been demonstrated experimentally by Dorn et al.[-l] (for a focused beam with 
no cavity). Parameters of the non-inset cavity and mode are R = 10, zi = 
0.05, and k = 7.2481 (A = 0.867). The cross section is at z = 0.3. The insets 
show the mode after it was followed to the perfectly hemispherical ^i = 
cavity shape; here the mode has an even smaller central spot and E^ ^ Ep. 
We expect a hemispherical mode to have only a single nonzero MB coefficient 
and this was verified by the solution: Oi = 1.0 and \ai\ < 1.1 x 10~^ for / 7^ 1. 
The hemispherical mode has k = 7.2243 (A = 0.870). 

ber of unknowns. When increasing /max or N^ij-s toward its problematic value, the A^ 
values for all k dramatically decrease, sharply lowering the contrast needed to locate 
the eigenvalues of k. The solutions that are found begin to attempt to set the field to 
zero in the entire cavity region, with some regions of layer outside the cavity having 
field intensities that are orders of magnitude larger than the field inside the cavity. In 
our experience this problem has not occurred in the two-basis method for /max near the 
semiclassical limit of/, Re{kR), and thus has never been a practical annoyance. On the 
other hand, the problem can occur at surprisingly low values of A^dirs- Yet taking iVdirs 
to be too low can often cause solutions to simply not be found: dips in the graph of 
Ar VS. Rek can simply disappear. The window of good values of A^dirs can be at least 
as narrow as A^dirs/10. The window of good /max values for the MB seems to be much 
wider (for many modes, it is sufficient to take /max to be half of Re{kR) or less). In this 
respect, the two-basis method is easier to use than the Bessel wave method. 

On the other hand, there were situations when a scan of A^ versus Rek with the 
Bessel wave method revealed a mode that was skipped in the same scan with the two- 
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Figure 11: A m = 2 mode for a conducting cavity {k = 7.885463). M2 is spherical and 
origin-centered with R = 10. Ml is a conductor; zi = 1.0. The inset shows 
ReEr^ in the x-z plane. 

basis method, due to the narrowness of the dip feature. To the best of our recollection, 
we have always been able to find a mode by both methods if we have made an effort to 
search for it. 

A direct comparison of the methods by looking at the solution plots rarely reveals field 
value differences greater than 1% of the maximum value when the modes are restricted 
to those that do not have a large high- angle component. The eigenvalues of k located 
by the two methods are usually quite close. Here we show one case in which there is a 
small but visible difference between the mode plots for the Bessel wave method and the 
two-basis method. 
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Figure 12: Top: Solution obtained with Bessel wave method, k 
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Bottom: Solution obtained with two-basis method, k = 8.57086 



zO.05900: A, 



2.9 X 10^5; / 



86. The field in layer X is not plotted. 



Figure 12 shows E^ in the x-z plane for a mode in a cavity with uq = 1.0 and nx = 0.5 
and no stack layers at all. This situation is meant to model a dielectric-filled dome cavity 
surrounded by air. (Our program assumes that layer is free space {uq = 1), so we have 
set nx < no to achieve this effect.) Here M2 is spherical and centered at the origin with 
i? = 10 and Ze = zi = 0.5. The reason for the disagreement between the two methods 
for this case is not known. This was the only low finesse cavity we have tried, as well 
as being the only cavity with a "high" index of refraction in layer 0. Another unusual 
property of this mode solution, which looks similar to a fundamental Gaussian, is that its 
electric field in the x-y plane actually spirals: its instantaneous linear direction rotates 
with z as well as with t. The other modes we show throughout the demonstrations 
section give extremely good agreement between the two primary methods. 

In our implementation, computation time to set up and solve Ay = b is of the same 
order for the Bessel wave method and the two-basis method (with variant 1). For the 
large cavity shown in Fig. 3, the Bessel wave method took about 10 s with A^dirs = 100 
and about 600 locations chosen on M2. Variant 1 with /max = 200, 600 M2 locations, 
and 600 uniformly spaced a^ , took about 15 s. For comparison, Variant 2 took about 
2 hours (achieving a 10~^ relative accuracy for each a^ integral). 
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5.5. Almost-Real MB Coefficients 

When Ml is a conductor we find that the imaginary parts of ai and bi are essentially 
zero. Furthermore, we often find that for dielectric Ml mirrors, the imaginary parts of 
ai and bi are one or more orders of magnitude smaller than their real parts. We show 
one way in which this tendency simplifies the interpretation of mode polarization in 
Appendix B. 

There is a an argument that suggests that ai and bi should be real for a conducting 
cavity. A conducting cavity has a real eigenvalues, k, and real reflection functions, 
T^s/p = — 1- It is not difficult to see that for real k the three types of M2 equations, 
E^ = 0, E\\ = 0, and H± = 0, can be satisfied with real ai and bi. The seed equation of 
course can also be chosen so that it can be satisfied by real variables. The Ml equations, 
because of the complex nature of fs/p = rs/pexp{—i2zikno cosa^ ) when zi ^ 0, are not 
obviously satisfied by real ai and bi. However, for a conducting cavity, the Ml boundary 
condition can be given as a real space boundary condition (like M2) instead of a Fourier 
space boundary condition. In this case, locations on Ml are chosen just as they are on 
M2. The same argument that the M2 equations can be satisfied by real ai and bi now 
applies. Thus the Ml boundary condition for a conducting cavity should not, in any 
basis, force ai and 6; to be complex. 

6. Conclusions 

The two methods presented in this communication have proved to be useful tools for 
modeling small optical dome-shaped cavities. As yet we have no general verdict on which 
method is best. Overall, it was quite helpful to have both methods available for finding 
all of the modes. Once a mode was found, agreement of both eigenvalue and eigenmode 
between the two methods was quite good. 

The primary advantage of both of these methods is the treatment of the dielectric 
stack through its transfer matrices. First of all, including the dielectric stack in the 
model is essential to obtain a number of qualitative features of the modes. A model 
with a mirror of constant phase shift is inadequate, except in special cases such as the 
fundamental Gaussian mode in its paraxial regime. In addition, including the stack in 
the model allows one to calculate the field inside the stack, which may be the location of 
interest for the engineer of an optical device. Secondly, the use of the transfer matrices is 
the most natural way to include a dielectric stack that is of large lateral extent compared 
to the wavelength. The only alternative that we know of would involve a separate basis 
for each stack layer, resulting in an increase of both dimensions of the matrix A by about 
A^ times, where A^ is the number of stack layers. 

Another advantage of the methods we have presented is their speed. Each point in 
Fig. 0, which was generated in a few hours or less, required Ay = b to be constructed 
and solved 20-60 times for a cavity with R ^ 50 A. The combination of using a basis 
expansion method, taking advantage of cylindrical symmetry, using one basis set for all 
the layers, avoiding numerical integrals, and using C++ inline functions for complex 
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arithmetic has resulted in a fast implementation. 

Several "stack-induced" effects that we have demonstrated here will be more thor- 
oughly treated in a separate publication. The splitting of the fundamental Gaussian 
(and other Gaussian modes) into a non-axisymmetric "V" shape may be of practical 
interest to workers who wish to achieve a tight focus. The m = modes, which have 
the greatest potential for tight focusing, are also of interest. The ability to analyze the 
persistent stack-induced mixing is a result that stands in its own right. This stack effect 
persists arbitrarily far into the paraxial limit; it exists as long as the resonance widths of 
both the cavity and the laser source are smaller than the splitting of the near-degenerate 
pair. This effect should exist for a wide range of cavity lengths, and is of interest to 
anyone engineering higher-order Gaussian modes. 

Perhaps the greatest limitation to the practical application of our methods is its rigid 
treatment of the curved mirror as a conductor. As mentioned in the Introduction, it is 
reasonable to expect highly focused modes, such as those shown in Figs. 3, 6 and 10, to 
undergo little change when the curved conducting mirror is exchanged for a dielectric 
one (with an appropriate change of mirror radius by —A/4 < 6 < A/4 to account for a 
phase shift). This is because the local wave fronts are primarily perpendicular to the 
curved mirror for such modes (imagine a Wigner function evaluated at the surface of 
M2). For paraxial geometries such as those in our discussion of stack-induced mixing, 
this argument says nothing about how replacing M2 with a curved dielectric stack will 
affect the modes. A pursuit of this question may require a more brute force approach. 

Despite all limitations, we find our methods to be extremely versatile and powerful 
(fast and allowing for relatively large cavities). The full vector electromagnetic field is 
used. Exactly degenerate modes can be separated. The shape of the curved mirror is 
arbitrary and we have implemented parabolic shapes as well as spherical shapes with 
no difficulties. A first-try scan of A^ versus Rek using either method with a non- 
contrived seed equation will locate the vast majority of modes that exist. Reasonable 
results are obtained when the interior of the cavity is a dielectric. Hopefully, the explicit 
development of the two methods given here can benefit a number of workers in optics- 
related fields. 

Note Added in Proof 

We note that we have not found significant stack-induced effects (the V mode of Section 
5.1 and the mixing of Section 5.2) for a "standard" quarter- wave dielectric stack which has 
design ABAB. . .AB with front surface A and ua > ub- Utia < nB, we do find both types 
of stack-induced effects. The standard {ua > tib) quarter-wave layer structure exhibits 
less variation in arg(rs/p(afc)) than our mirrors (cf. Fig. 15), and therefore behaves much 
like a perfect conductor except at high angles of incidence. The methods presented here, 
which model the curved dome mirror (M2) as a conductor, should therefore carry over 
to standard dielectric M2 mirrors. The planar mirror will in general have a "functional" 
layer design, leading to the stack-induced effects reported above. 
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A. Further Explanations and Limitations of the 
Model 

A.l. Exclusion of High-Angle Plane Waves 

As mentioned in the Introduction, in the usual application of a basis expansion method, 
the field is expanded in each dielectric region separately and henceforth each region gets 
its own complete basis and its own set of coefficients. However, our methods, when "cast" 
into the PWB, use a single set of plane waves, complete in layer 0, which are propagated 
down into the other layers via the transfer matrices. There are several advantages to this 
approach. The Bessel wave method becomes simpler, as there are far fewer unknowns. 
This approach is also ideal for incorporating the use of the MB for layer 0, as in the 
two-basis method, although a straightforward use of the MB for layer and a separate 
PWB for each layer g, g > 0, could be implemented. 

One drawback to our approach is that certain high-angle or evanescent plane waves 
are not included. No plane wave basis vectors are allowed for which UgSmO^ > riQ. If 
the true quasimode expansion in any layer has significant weight for these high-angle or 
evanescent waves, the calculated solution of the field could be erroneous. 

Probably the best way to determine whether or not such intrinsic error is present at 
a significant level is to look at the plane wave distribution in layer (using (58) or (82) 
to get this if one is using the two-basis method) and to see whether the distribution dies 
off as a^ approaches 7r/2. If it does, the solution should be reasonably error free. 

We note here that Berry [16] has shown how evanescent waves, in a finite region near 
the origin, can be expressed in the PWB (using plane waves with real-valued directions, 
as usual). Thus the fact that no evanescent waves are included in layer may not be a 
problem in itself for depicting the field in layer 0. 

A. 2. The Hat Brim 

A. 2.1. The infinitesimal hat brim 

For the vector field solution, a tiny hat brim must be added to M2 to correctly model 
a conducting edge. Figure 13 shows a cross section of the edge of M2 for increasingly 
better approximations of a real conducting mirror. In (a), Ep can be nonzero at the 
edge, which is unphysical. In (b), a small hat brim has been added and -E = at the 
inside edge, as it should be. In (c), locations are chosen on a (closed) conductor with 
finite thickness. As the thickness of M2 becomes insignificant compared to A, model (b) 
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Figure 13: Modeling the edge of the dome. 

should approximate model (c) arbitrarily closely. Usually 2-10 locations are chosen on 
the hat brim (the density of course is much greater on the hat brim than on the dome). 
We have taken Ub ~ O.OOOlyum for our demonstrations. 

A. 2. 2. The infinite hat brim and the 1-D half-plane cavity 

The way in which our model includes the interior of the cavity and the entire z > zi 
half-plane in a single region, layer 0, is somewhat unorthodox. We believe our method 
produces the correct field in a finite region surrounding the cavity, but not at z — > oo. 

A somewhat less strange model is obtained if one imagines an infinite hat brim. The 
upper half-plane is then no longer in the problem. Layer still extends infinitely in p 
as do the other layers, but the vertical confinement makes it easier to see that this is 
an eigenproblem and not a scattering problem. Unfortunately our implementation does 
not allow an infinite hat brim (unless Ml is a conductor and Ze = Zi, for which the hat 
brim condition is already set in the Ml equations). At least for some solutions, taking 
Wb ^ R and giving the hat brim roughly the same number of locations as the dome, 
produces no visible changes in the mode structure (from the solution obtained with a 
tiny hat brim). The mode shown in Fig. .3 is one of these. 

If we allow that some solutions are not affected by the value of Wb, we can take the 
infinite hat brim more seriously. The problem becomes similar in some respects to the 
1-D half-plane problem, in which x < is a conductor, x > L has a constant refractive 
index, and the [0, L] region is segmented into several 1-D dielectric layers. The solutions 
of the 1-D problem form a set of quasimodes, (or quasinormal modes), which are complete 
in [0, L] and obey an orthogonality condition on the same interval. The conditions for 
completeness, and the characterization of the incompleteness of the quasinormal mode 
for problems which do not meet the completeness conditions, are discussed in Ching et 
al. [17] and Leung et al. [18]. The model with the infinite hat brim does not rigorously 
meet the the completeness condition (for 3-D cavity resonators) . We note that it appears 
from the references mentioned above that the situation depicted in Fig. 14 would meet 
the completeness condition, with the quasinormal modes being complete in layers 0- 
A^. We do not attempt to resolve this further, as the use of our mode solutions as a 
basis for time-dependent problems is beyond the scope of what we have been trying to 
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Figure 14: Including the sides of the stack at radius Pfinite- 



accomplish. 



B. Negative m Modes and Sine and Cosine Modes 

The cylindrical symmetry group consists of 0-rotations and reflections about the x-z or 
y-z planes. For solutions with a fixed, nonzero m, rotations are equivalent to multiplying 
the field by a complex phase, which is of course equivalent to a translation in time. (This 
means that, when looking at any cross section plot in Section 5 for a fixed m 7^ 0, the time 
evolution can be realized by simply rotating the figure: counter-clockwise for m > and 
clockwise for m < 0.) Thus we say that for m 7^ 0, reflections generate a new solution but 
rotations do not. (If we were considering the sine and cosine modes discussed later in the 
section, instead of the ±m modes, we would say that a n/{2m) rotation generates a new 
solution but a reflection does not.) Hence the symmetry of the cavity causes modes to 
come in truly degenerate pairs. It can be shown that a reflection is equivalent to taking 
m to — m (up to a vr rotation). For m 7^ 0, general complex linear superpositions of the 
±m pair of modes are equivalent to arbitrary complex superpositions of any number of 
0-rotations, reflections, and combinations of these acting on a single +m mode. We note 
that since our methods solve y for a fixed m, we are able to separate the ±m pairs. 

The m = modes also come in degenerate pairs, but the interrelation can be more 
complicated than a refiection. The m = pairs can be separated easily in the two-basis 
method by choosing a seed equation that contains only ai or only bi coefficients, since 
the Ml and M2 equations do not couple ai and bi coefficients if m = 0. Thus all of the 
true (non-accidental), exact (not obtained only in a limit, such as the paraxial limit) 
degeneracies that exist can be separated. 

In constructing the Ml and M2 equations, it is easiest to write code that assumes 
that m > 0. For presentations and movies, it is nice to be able to plot both +m and 
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— m solutions, as well as the linear combinations adding and subtracting these modes. 
All four of these modes can be plotted from a solution y+m that has been been found 
using the positive azimuthal quantum number +m. Here we briefly discuss how to do 
this (for m 7^ 0). 

B.l. Plotting with -m 

To plot the — m modes, we use simple rules, allowed by the Ml and M2 equations, to 
create a new solution vector y^rn from the solution y+m- To obtain the held, the new 
coefficients from t/_m are simply inserted into the various expansion (or basis conversion) 
equations from Sections 3 and 4, using (— m) as the "m" argument in these equations. 
(All of these equations work for negative values of their m argument.) The relations 
Yi^.m = {-lT%*m, J-m = (-1)'"^™,, fi,-m = (-l)"^/im, and ^;,_„ = -(-l)'"^/,„ are 
useful. The rules for the new solution vectors are given in Table 2. (There is a certain 
arbitrariness to these rules, since (ay), with a being an arbitrary constant, satisfies the 
Ml and M2 equations.) 



Table 2: Rules for taking y+r, 



PWB 



MB 



y- 



scalar: 
vector: 



Su/d- 
Pujd,- 



^~1) Pu/d,+m 
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0'1,+m 
—bl,+m 



B.2. Plotting cosine and sine modes 

We can define the "cosine" and "sine" modes as X(c) = (X(+m) + (— l)'"X(_m))/2, and 
X(^s) = i^i+m) — (~l)™'^(-m))/2, where X stands for tp, E, or H. By adding explicit 
expansion expressions for the +m and —m modes, one can obtain expressions for X(c) 
and X(^s)- For example, using (61), (62), and the above rules one finds that 
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where ai and hi come from the original y+m- When a/ and bi are predominately real 
(see Section 5.5), we can see from the above equation that the real- valued physical fields 
Ep and E^ are proportional to cos m0 cos cut while E^ is proportional to sin m(/) cos cut. 
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(The cosine time-dependence of E is why we call this the cosine mode.) Thus for the 
cosine mode with \m\ = 1, the transverse portion of E has an average linear polarization 
along the x axis (note {Ey)^ = Vt), as opposed to the average circular polarization 
the m = ±1 modes would have. (For the m = +1 mode, {E • q)^ = Wt where 
q = — xsincjt + y cos out.) Due to their separating of time and (/'-dependence for the 
physical fields, the sine and cosine modes are very useful final forms of the field (when 
ai and 6/ are predominately real). 

C. Stacks used in Section 5 



Stacks I and II are similar to Ali_a;.Ga2,.As-AlAs stacks that Raymer has used experimen- 
tally [8]. Figure 15 shows the reflection phases for stack II. The stack parameters that 
are varied in our demonstrations are Ns and kg = 271 /Xg. The meaning of these param- 
eters can be inferred from the stack definitions below. The normal refiection coefficient 



for stack II with A^, = 20 is |r< 



/p[,ak 



0.9964. For A^, = 22, \r 



s/p\ 



0.9981. 
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Figure 15: Plane wave refiection phases (rad.) of stack II [Ns = 20) as a function of a^ 
(deg.). The solid (dashed) line is for s (p) polarization. The wavelength of 
the plane wave, Atest, is set at A^. The graph for stack I is similar. 



n 



{2N,) 



Stack I: no = nx = 1; rii = ris = . . . = n{2Ns-i) = 3.003; n2 = n^ = . . . 
3.51695. Layers l-(2A^s) are quarter- wave layers (optical thickness = As/4). 

Stack II: no = nx = I; rii = n3 = . . . = n(^2Ns+i) = 3.51695; n2 = ^4 = . . . = n(27v,) = 
3.003. Layers 2-{2Ns + 1) are quarter-wave layers. Layer 1 is a spacer layer that has 
optical thickness IXg. 
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